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ABSTRACT 

Angle-action coordinates are used to study the relic of an iV-body simulation of 
a self-gravitating satellite galaxy that was released on a short-period orbit within the 
disc of the Galaxy. Satellite stars that lie within 1.5 kpc of the Sun are confined to a 
grid of patches in action space. As the relic phase-mixes for longer, the patches become 
smaller and more numerous. These patches can be seen even when the angle-action 
coordinates of an erroneous Galactic potential are used, but using the wrong potential 
displaces them. Diagnostic quantities constructed from the angle coordinates both 
allow the true potential to be identified, and the relic to be dated. Hence when the 
full phase space coordinates of large numbers of solar-neighbourhood stars are known, 
it should be possible to identify members of particular relics from the distribution 
of stars in an approximate action space. This would then open up the possibility of 
determining the time since the relic was disrupted and gaining better knowledge of 
the Galactic potential. 

The availability of angle-action coordinates for arbitrary potentials is the key to 
these developments. The paper includes a brief introduction to the torus technique 
used to generate them. 
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1 INTRODUCTION 

Within the remarkably successful ACDM model, galaxy 
formation is a hierarchical process. Large galaxies, such 
as the Milky Way, are built up in me rgers and the 
accre ti on of smaller buildin g bloc ks (e.g. IWhite fc Reesl 
1 19781 : ISpringel fc HernquistJ |2003h . The signatures of 
these processes should still be visible today in the 
form of substructure such as streams in all compo- 
nents of the Milky Way (e.g. iLvnden-Bell fc Lvnden-Belll 



19951: iFreeman fc Bland-Hawthornll2002l : iHelmi et al.ll20o1 



Abadi et alj|2003h 



Evidence of the wealth of substructure in the stellar 
halo of the Milky Way has increased dramatically over the 
past 15 years, most notably wit h observations of the dis- 
rupting Sagittarius dwarf galaxy (llbata et al. 1994), and of 
the s treams visible in the SDSS data (e.g. Belokurov et al.l 
2006). Within the disc several substructures are known, as 
are several mechanisms that might be responsible for them. 
Stars are born in clusters within the disc, and over a pe- 
riod of several Galactic rotations these clusters evaporate 
and the stars become phase mixed, spreading in space but 
retaining closely related orbits. This is commonly referred to 
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as a supercluster - it is likely the Hyades- Pleiades superclus- 
ter formed this way l|Famaev et al.|[2005l ). Substructure can 
also be created by dynamical interaction with spiral arms, or 
a rotating bar component; for example the Hercules stream 
is thought to be associated wit h the bar of the Milky Way 
|Dehnenlll999l . l2000l : lFuxll200lf). I t has been suggested that 
the Arct urus group (lEggenl 1971! ) is debris from a merged 
satellite l|Navarro et al.ll2004r ). 

Many methods exist for finding this substructure, of- 
ten based on incomplete phase space information about the 
stars. In the outer parts of the halo it is possible to find 
substructure with only knowledge of stellar positions on the 
sky, sometimes in co njunction with photometric data (e.g. 
Belokurov et alj|2006|), or radial velocity measurements^ e.g. 



Lvnden-Bell fc Lvnden-Bellll995l : IHelmi fc Whitelll999D . In 



the solar neighbourhood, approaches that loo k for common 
prop er motions have been widely used, (e.g. IChereul et al.1 
Il999h . 

With the availability of full 6D phase space informa- 
tion for an increasing number of stars in the solar neigh- 
bourhood, most notably the catalogues r esulting from the 
Geneva-Copenhagen and RAVE surveys (iNordstrom et al.l 
120041 : ISteinmetz et all 120061 : IZwitter et al.ll2008l ), and with 
the prospect of a further increase by sever al orders of mag- 
nitude when Gaia data become available l|Perrvman et al.1 
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I2001I ), it is appropriate to consider methods for using these 
data in full. 

As discussed in iHelmi et all (|l999T ). the space of inte- 
grals of motion is a very promising one for finding substruc- 
ture such as superclusters or merger debris. Stars with a 
single small progenitor (star cluster or small galaxy) will 
have very similar values for the integrals, which will ensure 
they are tightly bunched in integral space even after phase 
mixing has produced a spatial distribution that is effectively 
featureless. There are additional benefits to using quantities 
that are not only integrals but adiabatic invariants, as these 
are more likely to remain constant as the Galactic potential 
changes over time. 

Previous work has focused on spaces defined by L z 
(which is an adiabatic invariant in an axisymmetric poten- 
tial), and other quantities that can be used as approximate 
integrals of the motion, such as the total angular momentum 
(only an integral of the mo tion in a spherically symmetric 
potential. IHelmi et al.lll999 | ). or are not adiab atic invariants, 
such as the energy ( Helmi fc de Zeeuwl 2000h or t he apocen- 
tre and pericentre of an orbit ( Helmi et al.1120061 ). 

In this paper we demonstrate the use of angle-action co- 
ordinates to find substructure in the solar neighbourhood. 
Actions are adiabatic invariants, and their conjugate vari- 
abl es, the "angle c oordinates", increase linearly in time. 
As lTremainej (l999) has pointed out, these properties make 
them exceptionally useful for analysing tidal streams. The 
difficulty of determining the actions of stars in non-spherical 
potentials has, however, severely restricted their application 
to date. This is the first in a series of papers in which we 
show how the concept of orb ital tori l|McGill fc Binne"vttl990l ; 
iKaasalainen fc Binneyfl994T l makes it possible to exploit the 
power of angle-action coordinates for practical galactic prob- 
lems. 

In Sections[2]&[3]we briefly introduce angle-action coor- 
dinates, and explain how we find them for stars with known 
phase space positions. In Section [3] we apply them to sim- 
ulated data of a satellite merger. We show that stars from 
the satellite that are observed in the solar neighbourhood 
are confined to a grid of patches in action space, even when 
rather a poor approximation to the Galactic potential is 
used. We show further that a diagnostic denned in terms 
of the angle coordinates enables the true potential to be 
distinguished from the false one. Section [5] discusses minor 
modifications to the analysis that are required to accommo- 
date secular evolution of the Galactic potential over the life 
of a relic. 



2 ANGLE-ACTION COORDINATES 

Three actions Ji and three conjugate angles coordinates 
8i provide canonical coordinates for s ix-dimensional phase 
space (e.g. iBinnev fc Tremaiiie] 120081 . §3.5). The conven- 
tional phase space coordinates w = (x, v) are 27r-periodic 
in the angles. The actions are conserved quantities for any 
orbit, and the angles increase linearly with time: 



0(t) = 0(0) + fi(J)i, 



(1) 



where the components of fi are the orbital frequencies. Thus, 
in six-dimensional (9, J) space, a bound orbit moves only in 
the three directions, over a surface that is topologically a 



three-dimensional torus. We generally refer to this model of 
the orbit as "the torus"; it is labelled by the actions J. 

Angle-action coordinates exist for any time- 
independent, integrable Hamiltonian. However, an analytic 
method of computing the transformations between normal 
phase space coordinates and angle-action coordinates 
w <-> (0,J) is only practical for the Hamiltonian de- 
fined by the isochrone potential, which as limiting cases 
includes the harmonic-o scillator and Kepler potentials 
l|Binnev fc Tremaindlioosl . §3.5.2). 

The Hamiltonian corresponding to a more realistic 
galaxy potential is not generally integrable. However, most 
orbits in an axisymmetric potential are approximately 'reg- 
ular' (non-chaotic), and thus admit three approximate iso- 
lating integrals of motion. Consequently, it is possible to find 
angle-action coordinates which describe motion on these or- 
bits over all interesting time-scales. 

Methods for constructing angle-action tori for orbits 
in a gener al potential have been in the literature for over 
a de cade l|McGill fc Binnevl Il990l ; IKaasalainen fc Binnevl 
1 1994T ) . but have been little utilised, primarily because of the 
technical challenges these methods present. It is, however, 
possible to encapsulate these technicalities so that users are 
protected from them, and once this has been done, it is 
nearly as easy to construct angle-action coordinates for an 
orbit in an axisymmetric potential as it is to numerically 
integrate the orbit with a Runge-Kutta routine, or similar. 

2.1 The torus method 

In our current implementation we restrict ourselves to orbits 
in axisymmetric potentials. Conservation of angular momen- 
tum, J^,, about the system's symmetry axis then reduces the 
problem to that of motion in the (R, z) meridional plane in 
the effective potential ^ e g(R, z) = &(R,z) + J|/2i? 2 (e.g. 
IBinnev fc Tremaind[200l . §3.2). J is the third action. 

We start with with a 'toy' Hamiltonian, H T , for which 
the relationship (R, z,pn,p z ) <-> (#|J, (9j, J^, jj) is known 
analytically^ namely that of a generalised effective isochrone 
potential 



*5i(r,0) = 



—GM 



+ 



b+ \fb 2 + (r - ro) 2 2 [(r - r ) sin & 



(2) 



where ■& is latitude (not to be confused with the dynamical 
angle coordinates); M, b, L z and ro are free parameters of 

the toy Hamiltonian. 

As described in detail in iMcGill fc Binnevl (|l99Ch . the 
toy torus J T = constant is distorted into the "target torus" 
that approximates the orbit by the generating function 



S(0 T ,J) =6 T ■ J + 2^S n (J)sin(n-0 T ), 



(3) 



where n is a two-vector with integer components and the 

1 In this paper we refer to variables such as the angles and actions 
in the toy Hamiltonian as 9 T , J T , and those in the target Hamil- 
tonian as 6 , J. This notation differs from that of lMcGill fc Binnevl 
lll990l) and IKaasalainen fc Binnevl lll994) . in which the toy angles 
and actions are referred to as 0, J, and those in the target poten- 
tial as 6', 3'. We make this change in notation as our focus is 
primarily on the application of this machinery, rather than on 
how it works. 
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notation n > indicates that the sum runs over exactly 
half the plane, with the origin excluded. The S n are free 
parameters of the generating function. The canonical trans- 
formation defined by this generating function is 



dS(0 T ,J) 







dS(0 T ,3) 
83 



J + 2^]nS n (J) 

n>0 



T ) 



= + 2 ^ dJ sin (n ■ 



1 



(4) 



(•>) 



(6) 



Since the transform is canonical for any values S n , and the 
toy torus is "null" in the sense that Poincare's integral J A dp- 
dq vanishes for any region A of the torus, the image of the 
toy torus under the canonical map (the target torus) is also 
null. Note that J is constant on the target torus, so in general 
J T is a non-trivial function of T . 

The values S n (J) (and the parameters of H T ) cor- 
responding to a given gravitational potential and J are 
found by enforcing the condition that the Hamiltonian H 
is constant on the target torus. Remarkably this condi- 
tion is sufficient to ensure that the torus corresponds to 
an orbit in the given galaxy potential. In practice we en- 
force this con dition by using th e Levenberg-Marquardt al- 
gorithm (e.g. IPress et al.l 1 1986T ) to minimize the statistic 
X 2 = jf- XX-ff — H) 2 , where the sum is over N p points 
spread evenly in T over the target torus. The derivatives 
of H with respect to the S n and the parameters of the toy 
Hamiltonian that the Levenberg-Marquardt algorithm re- 
quires can all be found through the chain rule. We generally 
refer to this process as an "action fit" . 

This minimisation determines the functional depen- 
dence w(# T , J). However, it does not tell us how w depends 
on 0: the S n have been determined for a single value of J, so 
dS n (3)/d3 is still undetermined. It is, however, possible to 
find approximate values for the frequencies, SI, by perform- 
ing an orbit integration over several periods in the target 
potential, starting from a phase space point on the torus, 
and performing a linear fit to the resulting values of Of (t). 

To find more accurate frequencies and expressions for 
the angle coordinates, we integrate small sections of the or- 
bit, startin g from a grid of points on th e torus, and use the 
equations ( Kaasalain en fc Binnev|[l994 ) 



0( O ) + nt = T (t) + 2j2^§y 1 



T (t)]. 



(7) 



Each integration for M time-steps yields 3M such equations 
in which 0(0), SI, and dS n (J)/dJ are unknowns. The equa- 
tions are linear in the unknowns and for M»l the number 
of available equations increases much faster than the number 
of unknowns. We truncate the sum over n to ensure that the 
number of unknowns is significantly less than the number of 
equations, and solve the equations using a least squares fit. 
We refer to this process as an "angle fit" . 



3 FITTING THE DATA 

The torus-fitting mechanism enables us to find w(0,J) for 
given values of J rather than enabling us to determine (0, J) 




Figure 1. Iterative procedure for determining J(w). We start by 
evaluating the planar orbit with energy E« (black dot). From its 
radial action and vertical frequency we move to the open circle 
(J R ,J' Z ). Further vertical moves are used to reach the line on 
which the energy is that of the star. Then we move along that 
line to the star, where the velocities of the orbit agree with those 
of the star. 



given the coordinates w of a star. We now explain how we go 
from w to (0, J) by an iterative procedure, which is some- 
what ad-hoc, but converges quickly. We have the numerical 
value L z of J^, so the problem is to find a point in the slice 
through action space — L z . Fig. Q] shows this slice. Sev- 
eral lines of constant H are shown: the full line is the line 
that corresponds to the energy E of the given star. The dot- 
ted line is for the "planar energy" _E" = |vfj + $ ef f(7?, 0) 
that the star would have if it were confined to the Galac- 
tic plane. Simple one-dimensional integrals enable us to find 
the action J R and time-averaged radius R of the orbit rep- 
resented by the intersection of the dotted line with the Jr 
axis. Our iterations start at this point, as J R is typically a 
good estimate for the true value of Jr. 

From there we move vertically up towards the full line 
by an amount J' z . Bearing in mind that Sl z = dH/d.J z , we 
estimate J' z from the first-order expression E ~ H(J R , J z ) ~ 
E n + Sl z J' z . That is, we take J' z = (E - E^)/v R , where we 
have approximated Sl z by the vertical epicycle frequency v 
at R. 

We obtain the torus ( J R , J' z , L z ) and using its energy E2 
and frequency Q z we obtain an improved approximation to 
Jz by incrementing J' z by AJ' Z — (E — Ezj/Slz. This step is 
repeated until the energy of the current torus is sufficiently 
close to E. 

Once we have converged onto the line H = E in Fig. [T] 
we can move along it with increments AJ to J that satisfy 
(to first order) CI ■ AJ = - in this process only the R and z 
components of J are changing, so A J is determined by AJ Z . 

If the orbit does not go through the location of the 
star, AJ 2 is increased when the maximum height at the 
star's radius is too small, and decreased when the orbit does 
not reach the star's radius. Once the orbit reaches the star, 
J z is adjusted until the local value of v z agrees with the 
observational value. 

This procedure has converged for all values of J that 
we have tried, (irrespective of, for example, whether J z is 
small) and typically involves ~ 20 torus fits per star. With 
our current torus-finding code (which we have not attempted 
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Figure 3. (left) and Q z (right) plotted against f2jj for a 

random sample of particles from our simulation. The lines shows 
the the rotational frequency of a circular orbit f2 c ; rc (left) and 
vertical frequency u (right) plotted against the epicycle frequency 



Figure 2. Particle positions in the x-y (Galactic) plane initially 
and after t = 1.5,4.5 and 7.5 Gyr (as labelled). 

to optimize for speed) the procedure requires ~ 15s per star, 
so even now actions could be found for tens of millions of 
stars in of order a week on a cluster of 1000 processors. 

After fitting a torus to the observations, we obtain the 
star's angle coordinates and more accurate frequency values 
by performing an angle fit. 



4 APPLICATION TO SIMULATION DATA 

To illustrate what can be achieved with the torus method, 
we use it to examine the debris of a self-gravitating satellite 
that was disrupted during an iV-body simulation in which 
the satellite moved in a fixed Galactic potential. We focus 
o n a case that is very similar to that described in Section 3.2 
of iHelrm" et all (|2006l ). which was designed to reproduce the 
properties of the Arcturus group. 

We represented the satellite by 5 x 10 5 particles in a 
King sphere of concentration c = log 10 (r t /r c ) = 1.25, core 
radius 0.39 kpc, and total mass 3.75 x 10 8 Mq. It was placed 
on an orbit in the plane of the Galaxy that has apocentre 
at 9.3 kpc, pericentre at 3.1 kpc, and angular momentum 
970kpckms _1 . The satellite is initially placed at apocentre 
and followed for ~ 9 Gyr. 

The sel f-gravity of t he satellite was found using 
GYRFALCON (|Dehnenl 120021 ) and the static Galaxy poten- 
tial was that of Model 2 in iDehnen fc Binnevl (|l998u . This 
model is axisymmetric, and consists of somewhat flattened 
spheroids representing the halo and bulge, and three expo- 
nential disc components to represent the gas disc and the 
thin and thick stellar discs. 

The particle positions in the x-y plane are plotted in 
Figure [2] for the initial conditions and at t — 1.5,4.5 and 
7.5 Gyr. After 1.5 Gyr the satellite is spread over all azimuths 
and particles are found over their entire radial range, but 
substructure is clearly visible. After 4.5 Gyr phase mixing 



has progressed to the extent that this structure is nearly 
undetectable from the physical positions alone. 

We calculated angles, actions and frequencies for satel- 
lite particles using both the true potential and a rather dif- 
ferent potential - a Miyamoto-Nagai potential (Section l4.3|) . 

4.1 The frequencies 

Fig. [3] shows a random sample of satellite particles at t — 
9 Gyr in the (fi^fiij) plane (left) and the (S} z ,Sln) plane 
(right) when the true potential is used. The full lines show 
the relationship between the epicycle frequency k ~ flu and 
the circular frequency f2 c irc — fi</> (left) or vertical frequency 
v ~ Q, z (right). These demonstrate that the strong corre- 
lation between and fi^> arises because each frequency 
depends strongly on energy and only much more weakly on 
either eccentricity or inclination to the plane. The relation- 
ship between £l z and the other frequencies is much broader, 
reflecting the strongly anharmonic nature of vertical oscilla- 
tions, which cause the vertical frequency to depend strongly 
on vertical amplitude. Because the simulated satellite is on 
an orbit in the Galactic plane, £l z is of little further interest 
in this study (though it will be in other cases). 

Consider now the frequencies of particles that lie in a 
given volume around the Sun, as survey stars usually do. 
We place the Sun 8 kpc from the Galactic centre along the 
line to the satellite's initial location and select particles that 
lie within 1.5 kpc of the Sun. Fig. [3] shows these stars at 
t — 1.5, 4.5 and 7.5 Gyr in the plane spanned by Q,pt/2n and 
Q.Rt/2-n - the number of rotations about the galactic centre 
and the number of radial periods, respectively. Now we see a 
clear substructure within the plot. Stars are found in patches 
at (close to) integer intervals in Q$t/2TT, and regularly in 
Q.Rt/2-n. 

The reason for this clumping is simple: to be in the so- 
lar neighbourhood at time t, having also all been near each 
other at an earlier point in time (when part of a single satel- 
lite), these particles must all have moved a certain amount 
in (j>, plus or minus an integer number of complete rotations 
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Figure 4. fl^t/27T (number of azimuthal periods) plotted against Qnt/2n (number of radial periods) for particles within 1.5 kpc of the 
"solar position" in our simulation after 1.5, 4.5 and 7.5 Gyr (left to right). The particles are separated into patches corresponding to those 
particles which have performed (approximately) an integer number of rotations about the Galactic centre, and are at the appropriate 
point in their radial oscillations. 



about the Galactic Centre. This is, in essence, a selection 
effect caused by taking a window of finite size. Even af- 
ter 7.5 Gyr, at which point phase mixing has rendered the 
spatial distribution essentially featureless (Figure [2}, there 
is manifest clumping in the f2 plot. In the case of Q.Rt/2-n 
the patches occur more frequently than at integer intervals 
because the orbits cross the radial range twice per radial pe- 
riod. The non-zero size of the patches reflects the non-zero 
size of the window, non-zero initial velocities of the particles 
relative to the satellite motion and the non-negligible mass 
of the satellite, which causes orbits to deviate from orbits in 
the Galactic potential at early times. There is also a small 
spread due to any errors in the value of f2 found - this is 
clearly a small effect as the patches are still distinct after 
7.5 Gyr. The spread in Q z among stars in these samples is 
no narrower than that of all the satellite's stars, and not sep- 
arated into patches, both because our window constrains z 
only weakly, and because the initial values of 9 Z range from 
zero to 2-7T as a result of the satellite starting from the plane 
with v z — 0. 



The number of patches in i"2-space increases approxi- 
mately as t 2 , since the number of integers that lie in the full 
range of £lRt/2ir or Q.^t/2-K is proportional to t (since the 
range of ft doesn't change). The size of individual patches 
is, to a first approximation, determined by the size of the 
window from which particles are chosen: if an orbit lies 
within the window over a range A9r (ignoring for simplicity 
the dependence on 9 Z , 9$), then in fi-space (as opposed to 
J7t/27r-space) each individual patch will have width A0R/t. 
Therefore the size of each patch in Jl-space is proportional 
to t~ 2 , so the total area of the patches is approximately 
time- independent. The other effects mentioned above that 
cause the patches to have finite size have a similar effect on 
the size of the patches with the exception of the error in 
measurements of Q. Some of the patches are restricted in 
size because they meet the edge of the envelope of available 
fi-values (the values found in the satellite as a whole - Fig- 
ure [3] left). This is probably most obvious in the left panel 
of Figure[4] in which one patch (at Q.^t/2-n — 8) is far larger 
than the others because it nowhere touches the envelope. 
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Figure 5. Jr plotted against J^> for a random sample of parti- 
cles at the beginning of the simulation (top) and after 9 Gyr of 
evolution (bottom). 
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Figure 6. Jr plotted against for particles within 1.5 kpc of the "solar position" in our simulation after 1.5, 4.5 and 7.5 Gyr (left to 
right). As in Figure [4] the particles are divided into patches which increase in number and decrease in size as t increase. 



4.2 The actions 



iHelmi et alj (2006) sought to identify substructure in the 
disc by calculating the locations of stars in "APL" space, 
which is the space spanned by apocentre, pericentre and 
J<£ (which they refer to as L z ). Apo- and pericentre can 
be considered to be integrals analogous to actions, so APL 
space is a mapping of action space. Hence it is of interest 
to examine the distribution of the satell ite's stars in action 
space for comparison with the results of IHelmi et al.l (|2006l ) 
although we shall find it less interesting than the frequency 
and angle spaces. 

Fig. [5] is a plot of Jr against for a random set of 
particles (reflecting the satellite as a whole) at the begin- 
ning of the simulation (top) and at t — 9 Gyr (bottom). The 
particles remain in the same general area of the (J^,, Jr) 
plane, but at early times the actions are not constant be- 
cause the satellite is self-gravitating and the strong negative 
correlation between Jr and seen in the initial conditions 
is replaced by a (rather weaker) positive correlation^ 

The strong negative correlation in the initial conditions 
arises because initially the satellite is at apocentre. A par- 
ticle that moves relative to the centre of the satellite in the 
opposite direction to the satellite's rotation about the Galac- 
tic centre has less angular momentum than one that moves 
in the opposite direction, and - in the absence of the satel- 
lite's self-gravity - would be on a more eccentric orbit, and 
thus have a higher value of Jr. 

The weaker correlation between Jr and J^ seen at 9 Gyr 
arises because the actions of a particle become constant 
when the particle is stripped from the satellite and starts to 
feel the latter's gravity only weakly. This occurs at pericen- 
tre, when the effect of combining motion within the satellite 
with the motion of the satellite is precisely opposite of what 
it is at apocentre. The extent of the correlation between Jr 
and J<t, at 9 Gyr is comparable to that shown in Fig. 5 of 
IHelmi et al.l l|2006l ) for the locations in APL space of stars 
that lie within 5 kpc of the Sun. 

Figure [6] is a plot of Jr against for particles that 
lie within 1.5 kpc of the Sun at t = 1.5,4.5 and 7.5 Gyr 



2 When the satellite's self-gravity is turned off, the actions prove 
to be constant as expected. 
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Figure 7. A projection in J space of the actions of particles 
within 1.5 kpc of the "solar position" in the simulation after 
7.5 Gyr (the same particles plotted in Figure [6] right). In this 
projection it is clearer that the stars are separated into individ- 
ual clumps in J space. 



(left to right). The actions of the particles found in the so- 
lar neighbourhood cover almost the entire range in J found 
in the satellite as a whole, but are separated into distinct 
patches that decrease in size and increase in number over 
time. The number of patches increases slightly faster than 
the area of each patch decreases with the result that by 
7. 5 Gyr the patches a re starting to merge into bands. Fig. 5 
of IHelmi et al.l l|2006l ) shows such a series of bands in APL 
space for stars that lie within 1.5 kpc of the Sun. 

Since ft is a smooth function of J, our study of the 
distribution of particles in frequency space explains their 
distribution in action space: the "allowed" values of J cor- 
respond to "allowed" values of ft, which are confined to 
patches. While the constraints on ft do not involve fl z , the 
constraints on J do depend on J z , because both Qr and fl^. 
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depend on J z . Consequently, the positions of stars in J-space 
form a relatively regular lattice, but the principal directions 
of that lattice are not parallel to the Ji axes. Therefore, the 
tendency of patches to run together in the extreme right- 
hand panel of Fig. [6] can be eliminated by plotting a differ- 
ent projection of action space. For example, Figure [7] is a 
plot of 0.507 Jr + 0.862 J z against J^> - a projection chosen 
by eye from a 3D visualisation of the distribution in J-space 
- and in this plot the patches are all distinct. In general, 
the patches will be most cleanly separated when the lattice 
is projected along one of it principal directions, since then 
points with the same Qr and Q$ but differing in il z are 
projected on top of one another. The optimum projection 
depends both on the potential and on the region of J-space 
occupied by the stars, but it can be straightforwardly iden- 
tified for any set of data because fl is found at the same 
time as J. 



MiyamoLo-Nagai potential 
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4.3 Working with incomplete knowledge 

In reality we do not know the Galaxy's potential a priori. 
In this subsection we show that satellite particles can be 
identified using even a poor approximation to the potential, 
and then the true potential identified from structure within 
the sample of satellite particles. 

We repeated the above analysis using orbital tori in the 
Hamiltonian for a Miyamoto-Nagai potential 



$mn(-R,z) = - 



GM 



^R 2 + (a + Vz 2 +b 2 ) 2 



(8) 



with mass M = 1.8 X 10 Mq, scale length a = 6kpc and 
scale height b = 0.3 kpc. This is a crude approximation to 
the true potential and one expects to be able to start from 
a better approximation to the Galaxy's potential. It is cho- 
sen such that the circular speed at the Solar radius is ap- 
proximately the same as in the true potential, and the scale 
height is similar to that of the true disc. We chose a scale 
length that is much greater than that of the true thin disc 
(which dominates the forces in the solar neighbourhood) as 
the Miyamoto-Nagai potential falls off quickly with radius, 
and we want to avoid any risk of having particles at or above 
the escape speed (at least one action diverges as a particle's 
speed tends to the escape speed). 

Fig. [5] shows plots of £Ir against Q$ (top) and Jr 
against (bottom). While the clear separation of particles 
into clumps in both ft and J seen in Figure [5] is somewhat 
smeared by using the wrong potential, it is not completely 
lost. Therefore, even when the true potential is unknown, 
these plots enable us to identify substructures. 

We can go further, and use the displacement of patches 
in f2 space by an erroneous potential to identify the true 
potential. Specifically, at time t the angle coordinates of the 
Qth particle 9 a satisfy 



fl a (t - t ) - (Oa - 9 a ,o) = 27rm a , 



(9) 



where the particle was at O a fi at time to; the (integer) com- 
ponents of m a give the number of oscillations that the par- 
ticle has made in R, z and <f>. 

There is negligible clumping in 9 Z because the satellite's 
orbit initially lay in the plane, but at some time to, before 
the satellite was stripped, individual values of 8R, a ,o were 
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Figure 8. Q^t/2ir (number of azimuthal periods) plotted against 
Q,Rt/2ir (top), and Jr against (bottom) for particles found 
within 1.5 kpc of the "solar position" in our simulation after 
4.5 Gyr as determined in the Miyamoto-Nagai potential described 
in Section 14.31 (for comparison see the middle panels of Fig- 
ures |1 & EJ . 



tightly correlated, as were values of 9$ :a ,o- In fact, in our 
simulation at t = 0, the satellite was centred at 8^^ — 
and at apocentre, where 9r ~ ir for all orbits. Therefore, we 
define the statistical measures 



i^R,at' — (0R ta — 9r,q) — 2-Km,R : a 



27rm<4,a 



(10) 



where the integers m^ a and mR tCt are chosen such that 5r iC , 
and are minimised; 6r : o = 7r and (J^o = (in this 
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Figure 9. 5r (solid line), and 8^ (dotted line) plotted against t' (see Equation 1 10 I t : the dashed line indicates the known true value of 
t. The lower panels (ii & iii) show 8r and 8^ determined using values of f2 and found in the same potential that the orbit integration 
was carried out in — the lower-right panel (iii) being a magnified section of the lower left hand panel (ii), focused around the true value 
of t. The upper panel (i) shows 8r and 8^, determined using values of fl and found in an incorrect potential (Equation [8}. There are 
strong minima in both 8r and 8a, around the true value of t when the true potential is used, whereas when an incorrect potential is used 
none is seen. 



case) and t' is a free parameter. In these expressions the 
numerators vary between ±7T, so when the values of fl a t' 
are randomly distributed, Si ~ 0.5, while when flit' has a 
well defined phase relative to the rest of the numerator, Si 
can approach either zero or 1, depending on whether the 
two halves of the numerator are in or out of phase with each 
other. 

At t — 4.5 Gyr in our simulation 5r and 8,/, were eval- 
uated by summing over particles within 1.5 kpc of the Sun. 
Fig. |9]shows the resulting plots of Sr and S^ as functions of 
t' with n a and 6 a determined in the Miyamoto-Nagai po- 
tential (upper panel) and in the true potential (lower panel). 
The upper panel shows that when the dynamical variables 
are evaluated in an erroneous potential, the Si scarcely move 
from their mean values after the first few megayears, and 
when one of them does move downwards, the other does 
not. By contrast in the lower panel for the true potential 
both Si display sustained beats as the relative phases of the 
two terms that make up the numerators in equations (|10|) 
have stable relative phases. In the neighbourhood of the true 



disruption time of the relic the beating swells in amplitude 
and, as the exploded view on the right of the figure shows, 
both move down together just 10 Myr past the relic's true 
age. 

When applying the test involving the Si to real obser- 
vational data, it will be necessary to search for beats over 
both t' and the mean of the initial azimuthal phases O^.c. 
Hence, Fig. [9] oversimplifies the problem because it shows a 
one-dimensional search rather than a two-dimensional one. 
Moreover, if the satellite was initially on an inclined orbit 
rather than one in the plane, the vertical angles 9 Z would 
be involved, so coincident beats would be required in three 
angle-dependent variables Si rather than two. Consequently, 
a smaller fraction of a relic's stars would satisfy this con- 
dition at a given time, making it harder to identify a relic 
against background noise. The upside of the involvement of 
8 Z is that it would give us the opportunity to constrain the 
vertical structure of the potential; a relic of a satellite that 
started from J z ~ only probes the Galactic rotation curve. 

Real data will also have background stars that did not 



Disassembling the Galaxy 9 



come from the satellite. Naturally, the Poisson noise in the 
distribution of background stars makes it harder to identify 
overdensities in J-space, but once the overdensities associ- 
ated with a remnant have been identified, background stars 
have negligible impact on the ability of the Si to determine 
the disruption time: any star in the overdensity is in the 
same part of f2-space as the remnant's stars, so the con- 
tributions to equation (|10[1 from background and remnant 
stars will differ only in the relevant values of 9. These differ- 
ences will inevitably be small because we are dealing with a 
small survey volume. For example, if 20% of the stars iden- 
tified as being part of the remnant are actually background 
stars, and are - on average - displaced from the position on 
the orbit where a typical remnant member with the same 
f2 would be by A8i — 0.2 (equivalent to ~ 1.5 kpc in the 
0-direction), this would only make a difference of 0.01 to 
Si. 

Before the satellite is completely disrupted, it is affected 
by dynamical friction, which will cause spreading in veloc- 
ity space. However this is a significantly smaller effect than 
that due to the self gravity and initial velocity dispersion in 
the satellite for an object of this mass - the Chandrasekhar 
dynamical friction formula (e.g. iBinnev fc Tremainel 120081 . 
§8.1) suggests that a satellite of this mass should be decel- 
erated by ~ 10kms~ Gyr - , and this decreases in direct 
proportion to the satellite mass as it is stripped. Hence in- 
cluding dynamical friction would shift individual velocities 
by significantly less than their intrinsic scatter, namely the 
internal velocity distribution of the satellite (18.6 kms -1 ). 

These matters, and any others that arise, will be inves- 
tigated further, but Fig. [5] clearly conveys the essential idea 
and gives a tantalising taste of the diagnostic potential of 
angle variables. 



5 SECULAR EVOLUTION 

Stars have continued to form at a significant rate throughout 
the lifetime of the Galaxy's thin disc, and it must be pre- 
sumed that the disc's mass has increased significantly over 
the last 5 Gyr. Since actions are adiabatic invariants, such 
secular evolution of the Galactic potential does not affect 
the distribution of relic stars in the — Jr plane (Fig. [5}. 
Secular evolution causes the frequencies at fixed J to be- 
come explicit functions of time, so we should write ft(J,t), 
and the increment in 8i over time t changes from A9i — flit 
to AOi = J^dtfli. The conditions for a relic star to be in 
the solar neighbourhood are given values of A9i mod 2n for 
i = R,(f>. Since with secular evolution Adi remains a contin- 
uous function of J, these conditions continue to be satisfied 
only in a grid of patches in action space; secular evolution 
shifts the patches, but does not blur them. Hence the di- 
agnostic power of plots like Fig. [6] is unaffected by secular 
evolution. 

In the presence of secular evolution it becomes necessary 
in equation (|9]) to replace f2 a (t— to) by J* dt tt a - To evaluate 
the required integrals, one must adopt a model of the history 
of the potential, which determines the time dependence of 
SI. The required model is self-evident if secular evolution 
is confined to growth in the disc's mass at a known rate. 
Uncertainties in this rate will make it harder to locate beats 
in Fig. H 



Of course the key to calculating the secular evolution 
of stellar systems, be they globular clusters or galaxies, is to 
express their distribut i on fu nctions in terms of actions (e.g. 
ISellwood fc McGaughl 120051 ). so the availability of orbital 
tori for arbitrary potentials opens up new horizons in this 
area. 



6 SUMMARY 

A major hope of "near-field cosmology" is to identify within 
the Galaxy groups of stars that were accreted together 
l|Freeman fc Bland-Hawthornll2002l ). We have demonstrated 
the power of angle-action coordinates for doing this by 
studying the debris of a self-gravitating satellite of mass 
3.75 x 10 s Mg, released within the plane of a realistic Galac- 
tic potential on an orbit with apocentre at 9 kpc. On this 
short-period orbit the satellite's stars become well phase 
mixed within a couple of gigayears, and they are quite widely 
distributed in action space. Nonetheless, the stars that lie 
within 1.5 kpc of the Sun are concentrated into a grid of 
patches in action space because only stars with certain fre- 
quencies are currently near the Sun. To see the patchiness 
of the distribution in action space it is not necessary to use 
the angle-action coordinates of the true Galactic potential. 
But the correct potential must be used if statistical mea- 
sures constructed from the angle coordinates of stars are 
to show a characteristic pattern of beats from which the 
time at which the relic was disrupted can be deduced. Hence 
our results suggest a two-stage procedure: first a reasonable 
approximation to the Galactic potential is used to identify 
relics through the clustering of their stars' points around 
the nodes of a grid in action space. Then once a relic has 
been identified, the Galactic potential would be adjusted un- 
til the angle-variable diagnostics showed pronounced beats. 
This second step would not only pin down the Galaxy's po- 
tential, but also reveal the time at which relic was disrupted. 

Growth in the mass of the disc since the satellite fell in 
would have significant effects only on Fig. U] to recover this 
plot it would be necessary to model the time dependence of 
the Galactic potential, so that the integrals J dtfl could be 
evaluated. We anticipate that with the help of angle-action 
coordinates this could be done to sufficient accuracy, but 
defer this refinement to a subsequent publication. 

We have neglected the deviations of the Galactic poten- 
tial from axis ymmetry. C o uld th ese deviations have a signif- 
icant impact? I Juric et~aT, I (|2008l ) use star counts in the SDSS 
survey to show that the Galaxy's thick disc is remarkably 
axisymmetric near the Sun. This finding suggests that it is 
legitimate to neglect the bar when searching for relics within 
the thick disc, such as the Arcturus group. In general, the 
quadrupole moments of the bar's gravitational potential will 
decline rapidly outside the end of the bar at R ~ 3 kpc, so 
stars that are not resonant with the bar will not be strongly 
affected by it. The observed axisymmetry of the thick disc 
suggests that few if any of its stars are resonantly trapped 
by the bar, so their orbits can be safely modelled with an 
axisymmetric potential. The question of how the - phase- 
dependent - effects of the bar would impact Fig. [9] may 
prove important. 

The exciting possibilities discussed here rest on two 
foundations. One is the availability of angle-action coordi- 



10 P. J. McMillan & J. J. Binney 



nates for any given potential, and the other is the availabil- 
ity of full phase space coordinates for significant samples of 
stars. The torus construc tion technique developed in a se- 
ries of paper starting with lMcGill fe Binnevl (|l99(t) can pro- 
vide angle-action coordinates, and programmes such as the 
Geneva-Copenhagen, RAVE and Gaia surveys will provide 
the phase space coordinates. 

Currently the torus technique is restricted to either ax- 
isymmetric systems or two-dimensional non-rotating bars. 
However, extension to three-dimensional bars, including 
bars that are rotating with a constant pattern speed, is in 
principle straightforward and will be attempted soon. 

Clearly when this technique is used to search a real cat- 
alogue for relics, and then to analyse them, difficulties will 
be encountered that we have ignored here. Most obviously 
one will have to contend with errors in the phase space co- 
ordinates of stars (primarily due to errors in distances) and 
with the difficulty in picking out overdensities in action space 
against a background of Poisson noise from field stars. We 
are currently applying the method to ~ 200 000 stars from 
the RAVE survey and hope to report the results in the near 
future. 
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